##库的导入
import pandas as pd
import numpy as np
from sklearn.decomposition import PCA
from scipy.stats import zscore
import scipy.cluster.hierarchy as sch
import pylab as plt
import seaborn as sns
from sklearn.cluster import KMeans
from sklearn.metrics import confusion_matrix
from sklearn.preprocessing import StandardScaler
from factor_analyzer import FactorAnalyzer as FA
from factor_analyzer.factor_analyzer import calculate_kmo
from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity
from sklearn.metrics import silhouette_score
import matplotlib.pyplot as plt
import libpysal
import pandas as pd
import numpy as np
import geopandas as gpd
import transbigdata as tbd
import esda
from pyproj import CRS
from splot.esda import moran_scatterplot
from sklearn.metrics import davies_bouldin_score
from sklearn import datasets
from sklearn.cluster import AgglomerativeClustering
###

#地理数据读取
gdp_data = pd.read_csv('D:/shuju/pingjia5.csv',encoding='ascii')  #导入数据
geography_data = pd.read_csv('D:/shuju/province_geography_data.csv',encoding='GB2312')
gdp_data2 = pd.read_csv('D:/shuju/pingjia0.csv',encoding='ascii')
data_stop = pd.merge(gdp_data, geography_data, on='province_name')
data_stop2 = pd.merge(gdp_data2, geography_data, on='province_name')
data_stop.head(30)
data_stop2.head(30)
###指标的处理
file_path='D:/shuju/shuju0.xlsx'
df66=pd.read_excel(file_path,sheet_name=5)
df11=pd.read_excel(file_path,sheet_name=0)
df22=pd.read_excel(file_path,sheet_name=1)
df33=pd.read_excel(file_path,sheet_name=2)
df44=pd.read_excel(file_path,sheet_name=3)
df55=pd.read_excel(file_path,sheet_name=4)
def defen(df):
    df=df.drop(['省/市地区','年份','乡村人口（万人）','耕地面积(公顷)','人均耕地面积（公顷）'],axis=1)
    # 按列取最大值、最小值,inplace=True)
    max = df.max(axis=0)
    min = df.min(axis=0)
    # 先全按照正向指标处理
    b = (df - min) / (max - min)
    # print(max.iloc[2])
    # 处理反向指标
    for i in [0, 2, 3, 5, 6, 7, 8, 9]:
        b.iloc[:, i] = (max.iloc[i] - df.iloc[:, i]) / (max.iloc[i] - min.iloc[i])
    # print(b)
    data = b + 0.01
    ###
    ####计算五个得分
    # 求出评价对象数目
    n = data.shape[0]
    # 求和
    s = data.sum(axis=0)
    # 求比重矩阵
    P = data / s
    # 求熵以及差异指数
    e = -(P * np.log(P)).sum(axis=0) / np.log(n)
    g = 1 - e
    # 求出权重
    w = g / sum(g)
    f = P @ w
    # print(P)
    # print(w)
    ind1 = np.argsort(-f)
    # print(ind1)
    ind11 = np.zeros(31)
    ind11[ind1] = np.arange(1, 32)
    # print(ind11)#输出排序
    indexx = ['产业兴旺', '生态文明', '乡风文明', '治理有效', '生活富裕']
    scound = P * w  # 第二指标得分
    chanye = scound.iloc[:, 0] + scound.iloc[:, 1]
    shengtai = scound.iloc[:, 2] + scound.iloc[:, 3]
    xiangfeng = scound.iloc[:, 6] + scound.iloc[:, 5]
    zhili = scound.iloc[:, 9] + scound.iloc[:, 7] + scound.iloc[:, 4]
    shenghuo = scound.iloc[:, 10] + scound.iloc[:, 12] + scound.iloc[:, 8] + scound.iloc[:, 11]
    combined = np.vstack((chanye, shengtai, xiangfeng, zhili, shenghuo))
    ai = pd.DataFrame(combined.T)
    d = ai * 1000
    #print(d)
    # 计算KMO与p值验证因子分析合理性
    kmo_all, kmo_model = calculate_kmo(d)

    ##bartlett球状检验
    chi_square_value, p_value = calculate_bartlett_sphericity(d)

    # 输出检验结果
    # print('卡方值:', chi_square_value)
    # print('p值:', p_value)

    r = np.corrcoef(ai.T)  # 相关系数矩阵（原指标）
    print(r)
    val, vec = np.linalg.eig(r)  # 求特征值
    cs = np.cumsum(val)  # 求特征值的累加值
    rate = val / cs[-1]
    srate = np.cumsum(rate)
    # print(srate)
    n = 3  # 主因子个数

    fa = FA(3, rotation='varimax')
    fa.fit(d)  # 求解方差最大的模型
    A = fa.loadings_  # 求解载荷矩阵
    gx = np.sum(A ** 2, axis=0)  # 计算方差贡献
    s2 = 1 - np.sum(A ** 2, axis=1)  ###计算特殊方差1
    ss = np.linalg.inv(np.diag(s2))
    f = ss @ A @ np.linalg.inv(A.T @ ss @ A)  # 计算因子得分函数的系数
    df = d @ f  # 计算因子得分
    pingjia = df @ gx / sum(gx)
    return pingjia,df,p_value,kmo_model

provinces=df66.iloc[:,0]
provinces=pd.DataFrame(provinces)
#计算各年份得分
pingjia0,yinzi0,p0,kmo0=defen(df11)
pingjia1,yinzi1,p1,kmo1=defen(df22)
pingjia2,yinzi2,p2,kmo2=defen(df33)
pingjia3,yinzi3,p3,kmo3=defen(df44)
pingjia4,yinzi4,p4,kmo4=defen(df55)
pingjia5,yinzi5,p5,kmo5=defen(df66)
data_arrays = [pingjia0, pingjia1, pingjia2, pingjia3, pingjia4, pingjia5]

# 合并数组
combined_data = np.column_stack(data_arrays)

# 计算平均值
mean_values = np.mean(combined_data, axis=1)

# 将平均值作为新的一列添加到合并后的数据中
combined_data=pd.DataFrame(combined_data)
combined_data['综合得分'] = mean_values
combined_data_sorted = combined_data.sort_values(by='综合得分', ascending=False)
combined_data_sorted['排名'] = combined_data_sorted['综合得分'].rank(method='first', ascending=False)
# 将numpy数组转换为pandas DataFrame，并添加列名
df_ = pd.DataFrame(combined_data)
df_.columns = ['2015年', '2016年', '2017年', '2019年', '2020年', '2021年', '综合得分']

# 基于平均值进行排序，并获取排名（注意：这里假设你想要降序排名）
provinces=df66.iloc[:,0]
provinces=pd.DataFrame(provinces)
# 将排名作为新的一列添加到DataFrame中
combined_data_sorted = pd.concat([provinces, combined_data_sorted], axis=1)
df_=combined_data_sorted
df_['排名'] = df_['排名'].astype(int)
df_.columns=['省/市地区','2015年','2016年','2017年','2019年','2020年','2021年','综合得分','排名']
# 打印结果
#print(df_)
df_.to_excel('paiming.xlsx')
pingjia5.to_excel('pingjia5.xlsx')
pingjia0.to_excel('pingjia0.xlsx')
###综合排名柱状图
df_sorted =df_.sort_values(by='综合得分', ascending=False)
#print(df_sorted)
plt.barh(df_sorted['省/市地区'], df_sorted['综合得分'],color='r')
# 设置标题和坐标轴标签
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.title('按序号排序的省份/市地区得分柱状图')
plt.xlabel('省/市地区')
plt.ylabel('综合得分')
plt.gca().invert_yaxis()
plt.show()
from scipy.cluster.hierarchy import linkage, dendrogram, fcluster
from scipy.spatial.distance import pdist
from sklearn import metrics
distance_matrix = pdist(yinzi5, 'euclidean')
Z = linkage(distance_matrix, method='ward')
plt.figure(figsize=(10, 7))
dendrogram(Z)
plt.title('Hierarchical Clustering Dendrogram')
plt.xlabel('Sample index')
plt.ylabel('Distance')
plt.show()
num_clusters = 4
labels = fcluster(Z, num_clusters, criterion='maxclust')
print("Cluster assignments:", labels)
silhouette_samples = metrics.silhouette_samples(yinzi5, labels)
print(silhouette_samples)



# 绘制样本的轮廓系数直方图
plt.hist(silhouette_samples, bins=50, range=[-0.1, 1])
plt.xlabel('Silhouette Coefficient Values')
plt.ylabel('Frequency')
plt.title('Silhouette Coefficient Histogram')
plt.show()
#计算每一类每项因子均值
N=yinzi5.shape[0]
yinzi5['category'] = labels
mean_values = yinzi5.groupby('category').mean()
print(mean_values)
mean_values.columns=['生活生态和谐因子','治理有效因子','乡风文明因子']
mean_values.to_excel('julei.xlsx')

# 分析明显有短板乡村
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# 绘制散点图
ax.scatter(yinzi5.iloc[:,0], yinzi5.iloc[:,1], yinzi5.iloc[:,2])
for i, txt in enumerate(yinzi5.index):
    x, y, z = yinzi5.iloc[i, 0], yinzi5.iloc[i, 1], yinzi5.iloc[i, 2]
    ax.text(x, y, z, txt, fontsize=5)  # 调整fontsize以适合您的图形
# 设置标签
ax.set_xlabel('生活生态和谐因子')
ax.set_ylabel('治理有效因子')
ax.set_zlabel('乡风文明因子')
plt.show()


#提取省份列
provinces=df66.iloc[:,0]
provinces=pd.DataFrame(provinces)
####乡村发展潜力
cp=yinzi5.max(axis=0)
cm=yinzi5.min(axis=0)
p1=np.linalg.norm(yinzi5-cp,axis=1)
p2=np.linalg.norm(yinzi5-cm,axis=1)
ff=p2/(p1+p2)
ff=pd.DataFrame(ff)
ff.insert(0, '省份/市地区', provinces)
ff.columns=['省份/市地区','得分']
ff_sorted = ff.sort_values(by='得分', ascending=False)
print(ff_sorted)
ff_sorted.to_excel('qianli.xlsx')



#空间
lat1sh=20.0347675
lat2sh=45.766287
lon1sh=87.5354511
lon2sh=126.6376322
bounds=[lon1sh,lat1sh,lon2sh,lat2sh] #设置地理区域的边界（经纬度）
params=tbd.area_to_params(bounds,accuracy=100)   #将地理区域边界转换为网格参数
#print(params)
data_stop['LONCOL'],data_stop['LATCOL']=tbd.GPS_to_grid(data_stop['lon'],data_stop['lat'],params)
data_count=data_stop.groupby(['LONCOL','LATCOL'])
#['lon'].count().rename('count').reset_index()  #按网格列和行索引对data_stop进行分组，并计算每个网格中的lon数量（可能是点的数量），然后重命名并重置索引。
data_count=data_stop.groupby(['LONCOL','LATCOL']).size().rename('count').reset_index()
data_count['geometry']=tbd.grid_to_polygon([data_count['LONCOL'],data_count['LATCOL']],params)  #将网格坐标转换
data_count['0']=data_stop['0']  #将data_stop中的'0'列复制到data_count中
data_count['province_name']=data_stop['province_name']
data_count=gpd.GeoDataFrame(data_count)  #将data_count转换为GeoDataFrame
data_count.head(40)

data_count.crs=CRS('EPSG:4326')
data_count_2416=data_count.to_crs(epsg='2416')  #将data_count的坐标参考系统转换为EPSG:2416
data=data_count_2416.copy()  #复制到data
k=10
wq=libpysal.weights.KNN.from_dataframe(data,k=k)  #使用libpysal库创建一个基于K最近邻（KNN）的空间权重对象，其中K设为10
#full_weights_matrix = wq.full()
#print(full_weights_matrix)
#print(wq)
#wq=libpysal.weights.KNN.from_dataframe(data,k=10)
#print(wq)
weight=data['0']  #将data中的'0'列赋值给weight变量。
lisa_moran=esda.moran.Moran(weight,wq)  #计算全局Moran's I统计量
lisa_moran.I,lisa_moran.z_norm,lisa_moran.p_sim  #获取全局Moran's I统计量的值、标准化Z值和模拟的p值
#print(lisa_moran.I)
lisa_moran_local=esda.moran.Moran_Local(weight,wq)  #计算局部莫兰指数
data['Is']=lisa_moran_local.Is
data['z']=lisa_moran_local.z_sim  #得到Is和z
data.insert(0, '省份/市地区', provinces)
data.to_excel('moran.xlsx')


# 绘制散点图
plt.figure(figsize=(10, 6))  #创建一个新的图形窗口，并设置其大小为10x6英寸
plt.scatter(data['Is'], data['z'], c=data['Is'], cmap='coolwarm', vmin=-1, vmax=1)  #绘制一个散点图
coeffs, cov_matrix = np.polyfit(data['Is'], data['z'], 1, cov=True)  #使用numpy的polyfit函数对data中的'Is'和'z'列进行线性拟合，并返回拟合系数和协方差矩阵。这里拟合的是一条直线
slope = coeffs[0]
intercept = coeffs[1] #从拟合系数中提取直线的斜率和截距
fit_x = np.linspace(min(data['Is']), max(data['Is']), 100)  #创建一个包含100个点的线性空间数组，这些点均匀地分布在data中'Is'列的最小值和最大值之间
fit_y = slope * fit_x + intercept  #使用前面得到的斜率和截距来计算线性拟合线上每个点的y值

plt.plot(fit_x, fit_y, 'black', label='Linear Fit: y = %.2fx + %.2f' % (slope, intercept)) #在散点图上绘制线性拟合线，颜色为黑色，并添加图例标签，显示线性拟合的方程
plt.colorbar(label='Is')  #显示颜色条
plt.title('Local Moran')
plt.xlabel('Local_Moran')
plt.ylabel('z')
for i, txt in enumerate(data_count['province_name']):
    plt.annotate(txt, (data['Is'][i], data['z'][i]))
plt.show()
plt.savefig("1.jpg")


###15年
lat1sh=20.0347675
lat2sh=45.766287
lon1sh=87.5354511
lon2sh=126.6376322
bounds=[lon1sh,lat1sh,lon2sh,lat2sh] #设置地理区域的边界（经纬度）
params=tbd.area_to_params(bounds,accuracy=100)   #将地理区域边界转换为网格参数
#print(params)
data_stop2['LONCOL'],data_stop2['LATCOL']=tbd.GPS_to_grid(data_stop2['lon'],data_stop2['lat'],params)
data_count=data_stop2.groupby(['LONCOL','LATCOL'])
#['lon'].count().rename('count').reset_index()  #按网格列和行索引对data_stop进行分组，并计算每个网格中的lon数量（可能是点的数量），然后重命名并重置索引。
data_count=data_stop2.groupby(['LONCOL','LATCOL']).size().rename('count').reset_index()
data_count['geometry']=tbd.grid_to_polygon([data_count['LONCOL'],data_count['LATCOL']],params)  #将网格坐标转换
data_count['0']=data_stop2['0']  #将data_stop中的'0'列复制到data_count中
data_count['province_name']=data_stop2['province_name']
data_count=gpd.GeoDataFrame(data_count)  #将data_count转换为GeoDataFrame
data_count.head(40)

data_count.crs=CRS('EPSG:4326')
data_count_2416=data_count.to_crs(epsg='2416')  #将data_count的坐标参考系统转换为EPSG:2416
data=data_count_2416.copy()  #复制到data
k=10
wq=libpysal.weights.KNN.from_dataframe(data,k=k)  #使用libpysal库创建一个基于K最近邻（KNN）的空间权重对象，其中K设为10
#full_weights_matrix = wq.full()
#print(full_weights_matrix)
#print(wq)
#wq=libpysal.weights.KNN.from_dataframe(data,k=10)
#print(wq)
weight=data['0']  #将data中的'0'列赋值给weight变量。
lisa_moran=esda.moran.Moran(weight,wq)  #计算全局Moran's I统计量
lisa_moran.I,lisa_moran.z_norm,lisa_moran.p_sim  #获取全局Moran's I统计量的值、标准化Z值和模拟的p值
#print(lisa_moran.I)
lisa_moran_local=esda.moran.Moran_Local(weight,wq)  #计算局部莫兰指数
data['Is']=lisa_moran_local.Is
data['z']=lisa_moran_local.z_sim  #得到Is和z
data.insert(0, '省份/市地区', provinces)
data.to_excel('moran15.xlsx')


# 绘制散点图
plt.figure(figsize=(10, 6))  #创建一个新的图形窗口，并设置其大小为10x6英寸
plt.scatter(data['Is'], data['z'], c=data['Is'], cmap='coolwarm', vmin=-1, vmax=1)  #绘制一个散点图
coeffs, cov_matrix = np.polyfit(data['Is'], data['z'], 1, cov=True)  #使用numpy的polyfit函数对data中的'Is'和'z'列进行线性拟合，并返回拟合系数和协方差矩阵。这里拟合的是一条直线
slope = coeffs[0]
intercept = coeffs[1] #从拟合系数中提取直线的斜率和截距
fit_x = np.linspace(min(data['Is']), max(data['Is']), 100)  #创建一个包含100个点的线性空间数组，这些点均匀地分布在data中'Is'列的最小值和最大值之间
fit_y = slope * fit_x + intercept  #使用前面得到的斜率和截距来计算线性拟合线上每个点的y值

plt.plot(fit_x, fit_y, 'black', label='Linear Fit: y = %.2fx + %.2f' % (slope, intercept)) #在散点图上绘制线性拟合线，颜色为黑色，并添加图例标签，显示线性拟合的方程
plt.colorbar(label='Is')  #显示颜色条
plt.title('Local Moran')
plt.xlabel('Local_Moran')
plt.ylabel('z')
for i, txt in enumerate(data_count['province_name']):
    plt.annotate(txt, (data['Is'][i], data['z'][i]))
plt.show()
plt.savefig("2.jpg")
